## IDTAXA for OTU table

args <- commandArgs(TRUE)

fasta <- args[1]
classifier <- args[2]
filename <- args[3]
threshold <- strtoi(args[4])
proc <- strtoi(args[5])
otu <- args[6]

require('DECIPHER')

# import the otu table
otu_data_frame <- read.table(otu, header = T, sep="\t")

## to rename the object as "trainingSet"
trainedClassif <- load(classifier)
trainingSet <- get(trainedClassif)

fastaString <- readDNAStringSet(fasta)
fastaString <- RemoveGaps(fastaString)

message("IDTAXA: performing the classification...")
ids <- IdTaxa(fastaString,
              trainingSet,
              type="extended",
              strand="top",
              threshold=threshold,
              processors=NULL, # any number seems to make R crash (macOS)
              verbose = F)

message("IDTAXA: classification done")

# extract assignments and confidence
assign <- array(NA, c(length(ids), 2))
colnames(assign) <- c("taxon", "idtaxa_confidence")
for (i in 1:length(ids))
{
  assign[i,"taxon"] <- paste0(ids[[i]]$taxon, collapse = ';')
  assign[i,"idtaxa_confidence"] <- ids[[i]]$confidence[length(ids[[i]]$confidence)]
}
# adding the seq ref and cleaning up
tmp <- c()
# if we used OTU filtering module
if (length(grep("cluster=", names(ids[1]))) > 0)
{
  if (length(grep("ASV", names(ids[1]))) > 0)  for (i in 1:length(names(ids))) tmp <- c(tmp, strsplit(names(ids[i]),split=";")[[1]][1])
  if (length(grep("ISU", names(ids[1]))) > 0)  for (i in 1:length(names(ids))) tmp <- c(tmp, strsplit(names(ids[i]),split="cluster=")[[1]][2])
} else {
  tmp <- names(ids)
}

rownames(assign) <- sub(";","",tmp) # if any (if OTUs generated by swarm or vsearch)
assign[,"taxon"] <- sub("Root;","",assign[,"taxon"])
assign[,"taxon"] <- sub("unclassified_Root","unassigned",assign[,"taxon"])
assign <- as.data.frame(assign)

## merging data frame with assignments, in the same order otu_data_frame[,1] works for both ASVs and OTUs...
output <- cbind(otu_data_frame,assign[as.character(otu_data_frame[,1]),])

write.table(output, file = filename, quote = F, sep="\t", row.names = F)
